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Abstract 



The freezing transition of hard spheres has been well described by vari- 
ous versions of density- functional theory (DFT). These theories should possess 
the close-packed crystal as a special limit, which represents an extreme testing 
ground for the quality of such liquid-state based theories. We therefore study 
the predictions of DFT for the structure and thermodynamics of the hard- 
sphere crystal in this limit. We examine the Ramakrishnan-Yussouff (RY) 
approximation and two variants of the fundamental-measure theory (FMT) 
developed by Rosenfeld and coworkers. We allow for general shapes of the den- 
sity peaks, going beyond the common Gaussian approximation. In all cases 
we find that, upon approaching close packing, the peak width vanishes pro- 
portionally to the free distance a between the particles and the free energy 
depends logarithmically on a. However, different peak shapes and next-to- 
leading contributions to the free energy result from the different approximate 
functionals. For the RY theory, within the Gaussian approximation, we es- 
tablish that the crystalline solutions form a closed loop with a stable and an 
unstable branch both connected to the close-packing point at a = 0, consis- 
tent with the absence of a liquid-solid spinodal. That version of FMT that has 
previously been applied to freezing, predicts asymptotically step-like density 
profiles confined to the cells of self-consistent cell theory. But a recently sug- 
gested improved version which employs tensor weighted densities yields wider 
and almost Gaussian peaks which are shown to be in very good agreement 
with computer simulations. 
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I. INTRODUCTION 



Some twenty years ago Alexander and McTague applied the formalism of Landau theory 
to the freezing transition of atomic materials 0. Using symmetry arguments they suggested 
that a bcc crystal should be the universally favored crystal structure, independent of in- 
teraction details. This theory attempts to describe the solid as a small, spatially periodic 
perturbation of a liquid. In a recent paper [|] we argued that such an approach should only 
be valid near the liquid-solid spinodal, at which the liquid state becomes locally unstable. 
The position of the spinodal is determined by the Fourier transform of the liquid direct 
correlation function c, and is given by the smallest density p for which the equation 

pc(p,k) = l (1) 

has a solution. Moreover, the perturbative approach does not apply to the local minima of 
the free energy in order-parameter space, which correspond to metastable or stable crystals, 
but rather to its saddle points. For the latter we confirmed universal behavior near the 
spinodal, which may have implications for nucleation 

The hard-sphere fluid has become the canonical model for freezing, since it captures in 
the most simple form the dominant packing effects while attractive interactions are believed 
to play only a secondary role. The best current theories for hard-sphere freezing are various 
versions of density- functional theory (DFT) 0-0. Usually they are explicitely constructed 
to reproduce the Percus-Yevick approximation cpy for the hard-sphere direct correlation 
function. In Fig. [I] we show the values of cpy(p,k) evaluated at the wave number k max (p) 
corresponding to the maximum at a given density p. One finds that there is no solution to 
Eq. (0) at physical densities p below the space filling density 6/7rx~ 3 where a is the particle 
diameter (at and beyond this limit cpy is not defined). This implies that those DFTs do 
not exhibit a liquid-solid spinodal at all. Therefore the saddle point solution branch of the 
stationarity equation derived from the density functional cannot connect to the liquid branch 
when the bulk density is increased. On the other hand, hard-core systems are characterized 
by a close-packing density as the maximum possible density of a given crystal structure. Upon 
approaching this limit a suitably defined crystalline order parameter, e.g., the inverse width 
of the density peaks, will diverge along the stable (minimum) branch. One may surmise 
that that is also true along the saddle point branch. Thus an alternative scenario to the 
bifurcation of a crystalline solution from the liquid at a spinodal point as discussed in 
Ref. 0, are two solid solution branches smoothly connected to each other at low densities 
which diverge at close packing and are completely isolated from the liquid. In order to test 
this hypothesis in the present work we examine the close-packing limit in detail using DFTs 
that have previously been applied to the low-density solid near the phase transition. 

Clearly, the strong localization of the particles in this limit provides an extreme case for 
such liquid-state based theories. Hence it is a good testing ground for assessing the qualities 
of different approximations. In contrast to most DFT studies of the hard-sphere solid we do 
not restrict the shape of the density peaks to Gaussians, but allow for general spherically 
symmetric peaks. This is especially interesting for the completely anharmonic hard-sphere 
crystal for which there is no a priori argument to justify Gaussians, even for small amplitude 
particle oscillations. 
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The starting point of density-functional theory is the free energy functional of the inho- 
mogeneous particle density p(r) with the general form 

F[p(r)} = F ld [p(r)} + FeM*)}- (2) 
The ideal gas contribution is given by [j3 = l/{k B T)} 

PFu{ P (r)] = f d 3 rp(r)[\np(r)X 3 - 1] (3) 

with the thermal de Broglie wavelength A. While the excess part F ex is not known exactly, a 
large number of approximate forms have been suggested and applied to various problems in 
the last decades [|]-[| . As we do not strive for completeness we will consider only two repre- 
sentative variants in this paper: the Ramakrishnan-Yussouff functional which is one of 
the first and simplest approximations that have been studied, and the fundamental measure 
functional developed by Rosenfeld and coworkers [fTQ|JTT|1 which at the present is believed to 
provide the best theoretical description of the hard-sphere fluid. From a given functional 
the equilibrium density distribution at a given bulk density pb is obtained by minimization 
under the constraint V~ l j d 3 p(r) = pb- The value of the functional at its minimum is the 
actual free energy of the system. For both functionals we performed numerical calculations 
at a series of bulk densities as well as an analytical analysis of the close-packing limit which 
enables us to determine the asymptotic density profile and free energy. 



II. RAMAKRISHNAN-YUSSOUFF THEORY 

A. Density functional and equilibrium profiles 

The Ramakrishnan-Yussouff functional follows from a density expansion of F ex around 
the homogeneous state truncated at the quadratic term: 

(3F ex /V = (3f ex ( Pb ) - ± J d 3 rd 3 r>(p(r) - p b )(p(r') - p b )c(p, |r - r'|). (4) 

Here f ex is the free energy density and c the direct correlation function (DCF) of the hard- 
sphere liquid at an effective density p, both of which are commonly approximated by the 
analytically known solutions of the Percus-Yevick integral equation. In a solid the density 
consists of a sum of identical peaks centered at the lattice sites R: 

p(r) = £> A (r-R). (5) 

R 

Throughout this paper it is assumed that the peaks are normalized 

J d 3 rp A (r) = 1 (6) 

and that the nearest-neighbour distance R nn in the lattice is determined by the bulk density, 
Rnn/c = (Pep/ Pb) 1 ^ 3 where a is the particle diameter and p cp is the maximum possible 
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density. In order to reduce the dimensionality of the integrations we moreover assume that 
Pa is spherically symmetric. Deviations from this symmetry exist 0,0, but are small 
especially near close packing fl4|. However, in contrast to most solid phase calculations 



which assume pa to be Gaussian here we do not restrict its shape. 
By insertion of Eq. (^) in Eq. (^) one obtains 

(3F ex /V = Pf ex (Pb) + \p 2 Ap, k = 0) - ^p b J2 [ drr 2 f dr'r' 2 p A (r)p A (r')w(r, r', R) 

r (7) 

where c is the Fourier transformed DCF and the integral kernel is given by 

/*2"7r /»1 /*1 

w(r, r', R) — 2tt / d(f) 12 / dcos9 / dcos9' 
Jo J-i J-i 

x c(p, (r 2 + r' 2 + R 2 + 2rRcos6- 2r'Rcos9' - 2rr'cos7) 1/2 ). (8) 

The angles 9, 9', and 7 are those between r and R, r' and R, and r and r', respectively, and 
cos 7 = cos 9 cos 9' + cos <pi2 sin 9 sin 9'. The contribution from R = simplifies to 

Hn 2 r r+r' 



w(r,r',0) = — -/ dr l2 r 12 c(p,r 12 ). (9) 
rr' J {r _ r ,\ 

Without loss of generality one may restrict the domain of pa to the Wigner-Seitz cell, so 
that the ideal contribution to the functional can be written as 

(3F t d/V = A7T Pb J drr 2 p A (r)[\np A (r)X 3 - 1]. (10) 

By minimizing and taking into account the normalization Eq. (^]) one finds the stationarity 
equation 

= exp[^ Eg / dr'r' 2 p A {r')w{r, r', R)] 
An J drr 2 exp[^ ^ R J dr'r' 2 p A (r')w(r,r', R)] 

The Percus-Yevick approximation for the hard-sphere DCF has the simple form 

c(p, r) = (c (p) + Cl (p)r + c 3 (p)r 3 ) Q(<r - r). (12) 

The density- dependence of the coefficients Cj can for example be found in Ref. \TE\. In the 
present context its most important feature is the cutoff at the particle diameter which leads to 
w(r, r', R) = for R — r — r' > a. Hence for the strongly peaked profiles in high density solids 
only the first shell of lattice vectors (|R| = R nn ) and the term with R = must be taken into 
account. We have calculated w(r, r', R nn ) by numerical integration using the trapezoidal rule 
with 50 3 mesh points, while an analytical expression for w(r, r', 0) was derived from Eq. @. 
The stationarity equation is then discretized in r and solved by iteration. An under relaxation 
scheme 

P (n+1) = wp&L + a - ")p {n) (13) 
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proved helpful to ensure convergence. Here the profile after the n-th iteration and 

pieL is the right hand side of Eq. ( |TTD calculated from p( n \ A typical value of the constant 
oo was 0.2. 

The resulting profiles are shown in Fig. |^. Their width scales with the free distance a = 
Rnn — <J that a sphere can move into the direction to its neighbour if the latter is kept fixed. 
The profile shapes approach a limiting form discussed below. Their most striking property is 
the occurrence of a maximum at intermediate distances r. This unphysical behavior vanishes 
in the close-packing limit. The DCF has been evaluated at the bulk density, p = p^. This most 
obvious choice has the disadvantage that the solid has a higher free energy than the liquid at 



all densities, as already pointed out in Ref. [16]. In the earliest DFT work the density of the 



coexisting liquid has been used instead, but that is not very reasonable when high density 
solids are considered. Other schemes to select a density p of an "effective liquid" have been 
proposed Ref. [fL6| , |r7] ,[3| , which always imply p < p\>. Figure [3] shows density profiles obtained 
with an arbitrarily chosen value p* = pa 3 = 0.95 which is close to the freezing density. 
Now the maximum does not occur and the convergence to the limiting shape is faster. The 
profiles are considerably flatter at small r than a Gaussian of the same width. 



B. Close-packing limit 

The results shown in Figs. ^ and |^ clearly demonstrate that, in spite of contrary claims 
TI||T9"|| , simple density-functional theories based on the Percus-Yevick DCF do exhibit a 



well-defined close-packing limit at which the peak width goes to zero. We will analyze this 
limit in more detail in the following. Let us assume that for small a = R nn — a the profile 

behaves as 

PA(r) = ^p (£) (14) 

with a width A = a/a where A, a — > with a fixed. We shall show that the stationarity 
equation has a solution consistent with these assumptions. The ideal free energy in this limit 
becomes (with N = p^V and s = r/A) 

POO 

(3 F id / N = An dss 2 p (s)[\np (s)-3\n(A/\)-l]. (15) 
Jo 

The relevant contributions to F ex are 

w(r = sA, r' = s'A, 0) = 16vr 2 c(p, 0) + 0(A) (16) 

and 

w(r = sA,r' = s'A,R nn ) = 2vr / d(f> 12 dx dx'c(p, a[l + A/a(a + sx - sV) + 0(A 2 )]) 

Jo J-i J-i 

= A7T 2 c(p,a)w{s,s',a) + O(A) (17) 

where 
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w(s,s',a) — — I ds 3 I ds' 3 Q(s' 3 — s 3 — a) 

ss J-s J-s' 

0, s' + s < a 

(s + s' — a) 2 /{2ss'), s' + s > a, —a < s' — s < a , , 

2(1 -a/s'), s'-s>a [ } 

2(1 — a/s), s' — s < —a 

Thus we finally have in leading order in A 

/*oo /»oo 

(3F ex /N = —2n 2 N nn c(p, a) / dss 2 / ds s' 2 po(s) po(s')w(s, s', a) + const = $ + const 

Jo Jo ( 19 ) 

where N nn denotes the number of nearest neighbours. 

The total free energy can now be minimized in two different ways. First, one can restrict 
to profiles of a fixed shape po(s), e.g., Gaussians, and differentiate only with respect to the 
scaled width a for fixed a which gives 

3 = -a^ (20) 

Due to the form of w(s, s', a) for a — > oo one has $ — > and thus the right hand side of 
Eq. (|^) also decays. On the other hand, for a — » $ tends to a positive constant (since 
c(p, cr) is negative), thus its derivative will be negative for sufficiently well behaved po( s )- 
Therefore the right hand side of Eq. (^Cj) is zero both at a = and a = oo and positive 
in between which implies a maximum at a finite value of a. This can be explicitely checked 
for Gaussians (po(s) = 7r~ 3 / 2 exp(— s 2 )) and step functions (po( s ) — 3/(47r)G(l — s)) for 
which the integrals in Eq. ( |l9l) yield |(1 — erf(a/y / 2)) and | — |a + |a 3 — -JjQ: 4 + ^a 5 . 



Depending on the height of this maximum Eq. ( p0|) has zero or two solutions. In the first 
case there are no stationary points with vanishing peak width at pb = p cp . This is the case 
for the "Onsager solid" discussed in Ref. which belongs to the same class of approximate 
functionals, but with c(p, r) replaced by its low-density limit — 0(cr — r). If — c(p, a) is larger 
(e.g., cpY(p[p C ,cr) = —20.345) the solution with smaller a corresponds to a saddle point 
and the solution with larger a to the stable solid minimum. We emphasize that the widths 
A = a/a for both solutions tend to zero for p& — > p cp . In Fig. [| we display the results 
obtained for fee and bec solids, employing Gaussian profiles and p = pb (fee: p* p = \/2, 
N nn = 12; bec: p* p = 3^/3/4, N nn = 8). We also include numerical solutions of dF/dA = 
for the nonasymptotic functional discussed above, evaluated for Gaussians. They approach 
the asymptotics quite slowly, especially for the saddle points. At low densities both branches 
are connected at an inflection point below which no solidlike solutions exist. 

Alternatively one can differentiate the asymptotic functional in Eqs. fll5"|) and (|19|) with 
respect to the profile po(s). Here one may set a = 1 without loss of generality. This leads to 
the Euler Lagrange equation 

p ^ = exp[irN nn c(p, a) / °° ds's' 2 p (s')w(s, s', 1)] 

1 ' 4n J °° ds s 2 exp[irN nn c(p, a) J °° ds's' 2 p Q (s')w(s, s', 1)] 

Its solutions, which represent the asymptotic profile shape, obviously only depend on the 
value of c at r = a, because near close packing the distance between two interacting particles 
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is always very close to a. The resulting shapes, shown in Figs. |^ and are rather flat close 
to the lattice site and decay strongly around r/a = 0.6, so they are definitely non-Gaussian. 

The iteration never converged to a second solution that would repesent the saddle point, 
even when started from the Gaussian saddle point discussed above. It has been conjectured 
in a DFT study of the isotropic-nematic transition of hard rods [^UJ that in general the 
saddle point is not accessible by iteration because it corresponds to an unstable fixed point 
(see also Ref. [gTJ). 

We mention a subtle point in connection with Eq. (|21|). Due to the form of w the right 
hand side goes to a constant for s ^> 1, which means that no normalized solution on [0, oo) 
can exist. However, as mentioned above, one may restrict to functions with a finite support 
(e.g., r < R nn /2, i.e., s < R nn /{2a)). For the numerical program indeed a much lower 
cutoff was used. In principle the solution now depends on the cutoff, but in practice this 
dependence is extremely weak because the constant approached for large s is of the order of 
exp(^N nn c(p,a)) ~ 10~ 53 so that the contributions from the tail of po(s) are neglegible for 
any reasonable value of the cutoff. Similar remarks apply to Eq. (|TT|). 

The free energy of the solid is determined by inserting the calculated equilibrium profiles 
into the density functional. Its asymptotic behaviour is given by 

(3F/N= -3\na + f + O(a). (22) 

The leading logarithmic contribution stems from and is in accordance with the result of 
free volume theory p2| and cell theory |]23|^^| . It has been proven exact for parallel hard 



cubes [^5j and for finite hard-sphere systems |26] and is generally believed to be exact also in 



the thermodynamic limit. The various theories differ in their prediction for the constant fo- 
In the Ramakrishnan-Yussouff approach (with p = p&) for an fee solid we obtain / = 21.7 
which is far above the molecular dynamics result fo = —1.493 [27|. As shown in Fig. |5| 



the asymptotic form is approached quite slowly, i.e., the higher order terms in Eq. 
are important up to high densities (which probably will also produce a bad equation of 
state). The free energies from the full minimization are only slightly below those for the best 
Gaussian profile (Fig. ]j|). 



III. FUNDAMENTAL-MEASURE THEORY 

A. Density functional 

Fundamental measure theory at present represents the best available DFT for strongly 
inhomogeneous hard-sphere fluids. In contrast to most previous approaches it does not de- 
pend on the direct correlation function as an input, but rather reproduces the Percus-Yevick 
correlation function as an output of the theory in the homogeneous limit. While the originial 



expressions [ 10] gave a divergent excess free energy for strongly localized particles, a recent 



empirical modification proved suitable also for the description of the freezing transition |TT 



We will call this version FMT1. Another new approximation has recently been suggested 



by Tarazona and Rosenfeld |28| based on more fundamental grounds. They presented a new 
derivation of FMT by enforcing the functional to reduce to exactly known expressions in the 
zero- and one-dimensional limit. They obtained a more complicated expression for one of 
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the excess free energy contributions that cannot be expressed in terms of weighted densities 
and also does not reduce to the Percus-Yevick free energy in the homogeneous limit. They 
also suggested a simplification by rescaling a certain expansion of this exact expression, 
which we adopt as FMT2. Due to its construction we expect FMT2 to provide a better 
description of the high-density crystal in which the individual particles are confined to quasi 
zero-dimensional cages formed by their neighbors. 

For a one-component hard-sphere fluid in three dimensions the fundamental-measure 
functional has the form 

3 

PF ex \p(r)]= d 3 rJ2Mn a (r)) (23) 

where the functions <pi depend only on the weighted densities 

n a {r) = J d 3 r'p(r)w a (r-r') (24) 

In FMT1 only two independent scalar and one vectorial weight functions occur: 

u>a(r) = 9(| -r), w 2 (r) = 5(| -r), w V2 (r) = -5(| - r), (25) 
for FMT2 a tensor weight function is necessary: 

*«(r) = ^*(f -r) (26) 
The expressions for the excess free energy density are: 

n 2 



ln(l-n 3 ), (27) 

2 ^2 



7TCT 2 



U 2 n V2 

27t<t(1 — n 3 ) ' 
[n\ - n 2 V2 f 



lFMTI _ 

3 " 24vm 3 (l -n 3 ) 2 ' 



(28) 
(29) 

iFMT2 _ 9 /on\ 

The density ansatz Eq. (|5J) induces a corresponding form for the weighted densities: 

^n^(r-R) (31) 



R 



with 



n { £\r) = J d 3 r'p A {r')w a {r - r'). (32) 

If pa is spherically symmetric the calculation of the weighted densities reduces to one- 
dimensional integrations: 
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7T 



r 
ira 



/r+a/2 / 2 \ / _ N ra/2—r 

dr'r' ^— - (r - r') 2 J p A (0 + 6 (- - rj 4vr y rfrV' 2 p A (r') (33) 

rr+a/2 

n A ( r ) = — / dr'r' p&(r') (34) 
ng^(r) = r +<T/2 dr'r' (r 2 - r' 2 + ^ p A (r'). (35) 

r r J|r— cr/2| V 4 / 

In this case the matrix n^(r) [defined by n(r) = n^(r- R)] is diagonal in any coordinate 
system aligned with r. An explicit calculation yields the eigenvalues 

rr+a/2 / ^2 \ 

n^^(r) = n A (r) = — r / dr'r' ( 4r 2 r' 2 — ( r' 2 — r 2 ) 2 ] p&{r') (36) 

2r y |cr/2— -r | V 4 / 

and 



nf > (r ) = ^ / dr'r' ( ^ - r' 2 + r 2 ) p A (/) . (37) 



-r+a/2 f ^2 

<\a/2- 

Note that Trn A (r) = n A (r). As p A is a strongly peaked function of width A the weighted 
densities n A (r), n A (r), and ^(r) have appreciable values only for |r — o~/2| < A while 
n A ( r ) tends to 1 for much smaller r and to for much larger r. Thus for small A at any 
point r in a solid at most two terms contribute appreciably to the sum in Eq. fl3~T|). 



We only consider fee solids. By exploiting the crystal symmetry the integration in Eq. (|23"D 
can be restricted to a simplex corresponding to 1/48 of the unit cell. In a coordinate system 
aligned with the conventional cubic unit cell its vertices are 

(0,0,0), ^(-=,0,0), i2 rw (— =,0), R 



y/T ' " ^" V 2V2'2V2' " ^" V 2V2' 2^2' 2V2 7 ' (38) 

It will be helpful to distinguish between the region A, that is "affected" by only one lattice 
site, and the region B affected by two sites, i.e. the set of those points whose distance to two 
sites differs from a/2 by a length of order A. As depicted in Fig. ^| region B consists of lens 
shaped sets around the midpoints between neighboring sites. Here the integrands 4>i(n a (r)) 
do not depend on the azimuthal angle around the line joining the sites, thus only a two- 
dimensional numerical integration over cylindrical coordinates p' and z' must be performed. 
In order to compute the full n in region B the contribution from one of the sites must be 
transformed to the coordinate system determined by the direction to the other site. This is 
accomplished by a rotation around an axis perpendicular to this direction by the angle 7 
given by 

co S7 = ^ = P'^^-KJ^ (39) 

r + r - [(/ 2 + (j' + fiW2) 2 )(p' 2 + (z'-fl„„/2) 2 )] 1/2 

Since in region A n a (r) depends only on the distance to the nearest lattice site the cor- 
responding integration can even be reduced to one dimension after the angular factors 
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stemming from the shape of the simplex have been worked out analytically. In practice 
a sufficiently large cutoff A (typically A ~ 2a) was chosen beyond which p A (r) is assumed 
to be zero, and the integrals over A and B were calculated separately. This approach proved 
to be much faster and more accurate than a straightforward 2d integration over the whole 
simplex, because then the integrand is essentially zero in large parts of the integration region. 



B. Equilibrium profiles 



In order to determine the equilibrium density profile under the constraint of spherical 



symmetry the functional derivatives of F ex are calculated. We first write 



5p A {r) 



d<pi 8n a (r') 
dn a 5p A (r) 



(40) 



and 



5n r 



5 



5p A (r) Sp A (r] 



R 



71 



(")/ 



fai">(d) 



R 



5p A (r) 



(41) 



d=r' R 



For n A '(d) the second term in Eq. (|33|) is rewritten as 0(er/2 — d)(l — J^/ 2 -d^ r ' r ' P&( r ')) 
which leads to 



5n 



(3) 



Sp A (r) 



a 



47rr i 9(- - d) 



©(M-|l 



(42) 



where we assumed that always r < a/2 + d. Furthermore one finds 



5n ( l\d) 
Sp A (r) 

5p A (r) 



ixar 



e(\d 



a 
V 



, 7rr 



a' 



d—[d 2 -r 2 + —)Q(\d-- 



(43) 
(44) 



For the tensor weighted density straightforward calculation leads to a similar but more 
lengthy expression. The partial derivatives d(pi/dn a are easily obtained from Eqs. fl27D — (p9|). 
The functional derivative can now be computed by inserting Eqs. (f42|)-(44) into Eq. (41) 
and that into Eq. (|40|) . For the integration over r' in Eq. fl40|) we adopt a similar scheme as 
for the functional itself. Due to the step functions in Eqs. (f42"D-(f44"D in region A the cutoff 



A can be replaced by the distance r for which the derivative is evaluated. In region B two 
terms from the lattice sum contribute. Because the integrand is nonanalytic at the lines 
where one of the distances d equals a/2 — r, a/2, or a/2 + r we partitioned the integration 
region B appropriately for the numerical integration. Together with the ideal free energy 
Eq. ( |i~0D one readily obtains the stationarity equation 



exp[ 



PA{r) 



1 5/3F ex /N i 
6p A (r) i 



4 7rr 2 



Ail J dr'r' 2 exp[- 



1 SpFex/N l 

'4?rr' 2 Sp A (r') i 



(45) 
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Again a mesh is introduced for p&{r) and the weighted densities are calculated by the trape- 
zoidal rule with linear interpolation between the mesh points. More sophisticated numerical 
integration routines are used for the integration over r' in Eq. (EDI) for regions A and B, and 
Eq. ( fi~5"D is iterated until the maximum relative change in p^{r) is less than 10~ 5 . 

The resulting profiles for FMT1 are displayed in Fig. [7| They are almost constant at 
small r and then decrease steeper around r = a/2, increasingly fast upon approaching the 
close-packing limit. In the next section we show that the limiting shape indeed is a simple 
step function. The profiles for FMT2 shown in Fig. |8| exhibit a much smoother, Gaussian-like 
decay and their width, measured, e.g., by (r 2 ) = J d 3 r r 2 p&(r), on the scale a is considerably 
larger than for both the RY and FMT1 functionals. Clearly, again the absolute width goes 
to zero linearly with a, as expected in the close-packing limit. 

C. Close-packing limit 

As for the Ramakrishnan-Yussouff functional we assume that asymptotically the density 
profile has the form Eq. flU] ) with A = a = aS. We have seen that in the important range 
the argument of the weighted densities is close to a/2. Therefore we set r/a = 1/2 + tS and 
determine the leading contributions to n^\r) for small 5 and fixed t: 



t=0 

"OO 



( 2 ), 

a 

i=0 



2tt ' x 



26(-t) / dss 2 p (s) 
Jo 

+ 1 dsspo(s) ((s-t) + 5(t 2 -s 2 ) + 25%s 2 -t 2 ) + ---) (46) 
J\t\ 

1 oo 

) (t) = -J2n2>(t)5*- 1 

i=0 
poo 

/ dssp (s) (l-2t5 + 4t 2 5 2 + ■■■) (47) 
J\t\ 

n { P ] (r, t) = r^- / dssp (s) (l-2tS+ (6t 2 - 2s 2 )6 2 + • • ■ ) (48) 
<™ J\t\ 

POO 

n W ( t ) = AnS / dsspo(s) (s 2 -t 2 + ---) (49) 
J\t\ 

n poo 

nf ] (t) = =f / dssp (s) (1 - 2td + 4(2t 2 - s 2 )5 2 + ■ ■ ■ ) (50) 
<5 J\t\ 

where the caret denotes a unit a vector. Since for any po(s) the first two terms in the 
expansions of and |n^ 2 ' ) | are identical one has n 2 — riy 2 = 0(5°) in region A. On the 
other hand, in region B the contributions to ny 2 from the two lattice sites have almost 
opposite directions so that n 2 , — riy 2 ~ 5~ 2 there. For FMT2 we find det n ~ 5 in both 
regions, because, due to the quadratic dependence of w»j(r) on the components of r, the two 
contributions do not cancel each other in region B. Taking into account that the volumes 
of A and B are proportional to 5 and S 2 , respectively, we can estimate the order of the 
individual free energy contributions $j = iV -1 J d 3 r<j)i(n a (r)): 
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A B 
$1 5° 5 
$ 2 S 5° 
$™ T1 5 A 5- 1 
$™ T2 5 2 5 3 

Thus at this point a qualitative difference between the two approximations arises, as 
different terms become dominant in the close-packing limit. We first discuss FMT1, for 
which $3# is the leading term. In a cylindrical coordinate system (z' , p' , <ft') , centered at the 
midpoint between two sites and with its axis directed towards (see Fig. |5|) one of them, the 
distances r± to the sites, which occur as the argument of the weighted densities are 



.'2 //r 2 



p' 2 + (z'±R nn /2) 2 } 1/2 . (51) 



In the scaled coordinates p = p' /(5o ) and z = z' /(5a) one has t± = 1/2 + p±z + 0(5) and 



4 

n l - n 2 V2 = ^n 2 o(t + )n 2 o(t-) H (52) 



which finally yields 



gfMTi „ 32 f°° f n 20 (t + )n 20 (t_) \ 1 

3B ~ 5 J dp J dZ \n w (U)+n 20 (t„)) (1 - n 30 (t + ) - n 30 (t_)y [M) 

Since n 3 o(t) G [0, 1] and n 2 o(t) > this expression is positive. It attains its minimum value 
zero for all profiles po(s) that have a strict cutoff at s = 1/2 so that region B is empty. In 
this restricted class of profiles the dominant contributions are $ia and F^. The former can 
be written as 



yl/2 

$ia = - dtn 20 (t) ln(l - 7i3o(*)) + 0(5). (54) 

J -1/2 



But the fact that n 2 o(t) = —dn 30 /dt implies $1,4 = 1 + 0(5) for all profiles. Since here the 
peaks around different sites are independent of each other, this result is consistent with the 
extensively discussed 0D limit of the fundamental-measure functional [11,^8|: For density 



profiles /?a(f) constrained to a volume that cannot hold more than one particle the exact 
excess free energy is (3F ex = 1 if J d z rp&(v) = 1. One of the merits of the present theory is 
that this limit is almost exactly fulfilled [|11|]. At last we are left with the ideal free energy 
Eq. (|i~5l) as the only relevant 0(5°) term, which, naturally, favors an evenly distributed 
density: 

Po(«)--e(i- s ). (55) 

7T 2 

This finding implies that the usually assumed Gaussian peaks represent a particularly 
bad approximation in this case. Indeed, in the Appendix we show that the width A of the 
best Gaussian is asymptotically related to the free distance a by a ~ A^y\n(— A) which 
means that the ratio A/a tends to zero, albeit very slowly. The intuitive reason is that the 
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tail of the Gaussian profile leads to an unfavorable free energy contribution $35 that can 
only be kept small if the tail increasingly "retracts" . 

Actually the above arguments for the asymptotic step function shape in FMT1 can be 
generalized to nonspherical profiles. Starting from 



PaW = -^p (r/A) 
and setting again s = r/A, A = a5 = a, and |r|/<r = 1/2 + t5 one has 
n { l\r,t) = J d 3 sQ(rs-t-5(t 2 - 2trs + s 2 ))p (s) 
Expanding for small 5 gives 

n { l\r,t) = J d 3 sQ(rs-t)p (s) + O(5). 

Analogously we find 



(56) 



(57) 



(58) 



r4 2) (r,t) = / d i sp (s 



5 



5(rs -t)- (t 2 - 2trs + s 2 )5'(rs - t) + 0(5) 



(59) 



and, using n^ 2 \r) = -Vn^'(r 



(3), 



n ( l\r,t) = I d 6 sp (s 



5 



r5(is - t) - r (if - 2trs + s 2 )5'(rs - t) 



-2(s- (fs)f)5(fs - t) +0(5) 



(60) 



Since the last term in this equation is perpendicular to f the combination n 2 , — riy 2 is 
still of order 5° in region A. If in region B the same coordinates (z, p, <$') as before are 
used and the vectors to the nearest lattice sites are denoted by r±, the fact that r + r_ = 
-1 + 0(5) yields 

n 2 ~ n V2 = 0(5 2 ). Thus, in summary all estimates for the individual 
terms given in the table above remain valid. Again the dominant term is positive and 
minimized by cutoff profiles. As the leading terms of the scalar weighted densities are related 
by ^2o(r, t) = —dn 30 (r, t)/dt the contribution $ ly 4 in leading order is still independent of the 
profile. The ideal term now enforces p&( r ) to be constant in the maximum allowed region 
C that is compatible with B = 0. It can be constructed by shifting the bounding planes of 
the Wigner-Seitz cell inward by cr/2 (see Fig. A given point r in C contributes to the 
weighted densities at r' only if |r — r'| < a/2. By construction all such r' lie within the 
Wigner-Seitz cell and thus cannot be "reached" from any r in the cell C around another 
site, which means that B is indeed empty. However if a point P outside of C were added 
the distance to its mirror point P' with respect to the closest Wigner-Seitz boundary plane 
would be less than a so that elements of B would lie on their joining line (see Fig. |S|). The 
cell C constructed here is identical to that of the self-consistent cell theory |23|. Its volume 
for an fee solid is a 3 / \/2. 
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We now turn to the second approximation (FMT2) for which $2B and are the domi- 
nant contributions. (Remember that $1^ is independent of po(s) i n leading order.) Analogous 
to Eq. ( |53"D we have, up to higher orders in 5, 



00 poo 



2B - 24 / dp dz- — — - 61 

1 - n 30 (t + ) - 7130 (*_) 



'0 JO 

with t± = 1/2 + p ± z. The corresponding stationarity equation is 

1 S&2B 



P V i7TS 5 po(s) 7 

p (s) = 7 r-. (62) 

47T j Q as s exp ^ i^ 5po(s /) J 

The functional derivative is calculated as in Sect. |111 BL The resulting asymptotic profile 



shown in Fig. |8] is close to those obtained for finite densities using the full functional. It is 
almost, but not exactly Gaussian. 

The asymptotic free energy of the fundamental- measure theory also has the form Eq. (p2|). 
In FMT1 one has _/ = ^ In 2 = 0.3466. Under the constraint of spherical symmetry this is 
replaced by fo = ln(6/7r) = 0.6470, both of which are much closer to the correct value 
than the Ramakrishnan-Yussouff theory. Results for finite distance from close packing are 
given in Fig. [| and agree, probably by accident, rather well with the computer simulations. 
The approach to the asymptotic law is quite slow. On the other hand for FMT2 not only 
the profiles but also the free energies (Fig. [5]) approach their asymptotic limit faster in this 
version of the theory. The value of the constant in Eq. ([22]) is found to be fo = —1.527 in 
very good agreement with the MD results. However, in view of the relatively large change 
in fo due to nonsphericity of the profiles as observed for FMT1, this may well be fortuitous. 
We did not consider nonspherical profiles in FMT2. 



D. Saddle point 



In view of the discussion in the introduction it would be interesting also to keep track of 
the saddle point between the liquid and the solid state when close packing is approached. 
Unfortunately, again one is plagued by the fact that the iteration of the stationarity equa- 
tion does not converge to a second solution. Furthermore the arguments of the asymptotic 
analysis do not apply to the saddle point because they essentially involve a minimization in 
two steps. Hence one must revert to parametrizations of the density with a few parameters. 
For Gaussians and step functions the saddle point occurs at a width A proportional to a 2//3 
in FMT1 and to a 1//2 in FMT2. However, a priori there is no reason to assume that at the 
saddle point the profile has a similar shape as at the minimum. We also tried profiles of 
the forms p&(r) ~ exp(— (r/A) n ) and p&(r) ~ (1 + r/A) _n and found numerically that in 
both cases the free energy at the maximum with respect to A decreases with decreasing n, 
down to the lowest feasible values of n. This suggests that the actual saddle point profile 
may decay very slowly, while within these restricted classes of profiles a true saddle point at 
a nondegenerate profile seems not to exist. 
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IV. MONTE CARLO SIMULATIONS 



Although an extensive computer simulation study of the density distribution in hard- 
sphere crystals has been carried out before |LJJ, no useful results for the radial distribution 



function have been published. In order to assess the quality of the various theories we 
therefore undertook a small Monte Carlo (MC) simulation ourselves. In an NVT ensemble 
of 8 3 spheres in an fee arrangement we measured the distribution of the particles' distance 
r from their equilibrium sites. We corrected for the movement of these sites due to shifts 
in the center of mass. Measurements were taken over 2 • 10 6 MC steps per particle for 
two bulk densities. The results are plotted in Fig. [10] on a logarithmic scale versus (r/a) 2 
and compared to the various DFT calculations. The quantitative agreement is excellent for 
FMT2. The profiles are close to Gaussians but decay faster at large distance than a Gaussian 
fitted to the small distance part. The dependence of the scaled profiles on bulk density is 
rather small in the examined range, but still qualitatively reproduced by the theory. The 
actual width of the profiles will increase with increasing particle number [I4j , but we did not 



attempt to correct for finite size effects. In Ref. [|T4| it was found that in the thermodynamic 
limit for high densities the width behaves as (r 2 ) 1 ^ 2 /a = 1.098 ± 0.004, again in almost 
perfect agreement with FMT2, for which (r 2 ) 1 ^ 2 /a = 1.025. This means that FMT2 is the 
first DFT which yields the correct value of the Lindemann parameter. 



V. SUMMARY AND DISCUSSION 

In summary, we have analyzed the close-packing limit of the hard-sphere crystal using 
three versions of DFT. All of them predict a peak width A that vanishes proportional to the 
free distance a and yields a logarihmic term in the free energy (see Eq. fl2"2p) stemming from 
the ideal gas entropy. Numerically this has been observed before, for Gaussian peaks, in 
two other DFTs, the generalized liquid approximation (GELA) and the modified weighted 
density approximation (MWDA) ||29|| . For the latter, however, it was found later that the 
solutions correspond to "unphysical" branches |3(| . The relative performance of the different 



theories can be judged from the profile shape obtained by free minimization. RY gives too 
narrow profiles with an unphysical maximum if the bulk density is used as the expansion 
point (Fig. 0). The shape and width are also wrong for other expansion points (Fig. 
FMT1 predicts asymptotically steplike profiles confined to the cells of cell theory (Fig. 
Only the FMT2 profiles (Fig. |8|) are in quantitative agreement with simulations at high 
densities (Fig. [10]) . In spite of the anharmonicity of the hard-sphere crystal they are close 
to Gaussians. Similarly the results for the next-to-leading free energy contribution improve 
from RY to FMT1 to FMT2 (Fig. |5[), the two FMT versions being much closer to the 
correct result than RY. This could have been expected from the way the RY approximation 
is constructed: A density expansion around a liquid state certainly is difficult to justify for 
the highly ordered high-density crystal. 

If one restricts the profiles to a fixed shape saddle points of the free energy are found 
at widths decaying ~ a x with xry = 1, xfmti = 2/3, and xfmt2 = 1/2. Insofar the 
global scenario for the crystalline solutions proposed in the introduction is comfirmed (see 
also Fig. |j). However, as detailed in Sec. [Ill D| , in larger classes of functions the saddle 
point remains elusive. We remark that the saddle point is a property closely connected to 
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the mean-field type free energy functional and, e.g., is not directly accessible by computer 
simulations. 

Comparing the two variants of FMT we see that the structure in the close-packing 
limit is sensitive to subtle differences between DFT approximations and thus might be a 
guiding line in the construction of better FMT-like functionals. Besides our FMT1 some 
other approximations for 



have been suggested in Ref. 



with 



n[ which are all of the form 
ny 2 



n 2 



(63) 



T ° (l-n 3 ) 2 ' 

The power of n 2 is determined by dimensional arguments and the function / can only 
depend on the absolute value of £ because of the isotropy of space. From Eqs. fl47f ) and 
( fPf) we have £ = 1 + 0(5 2 ) in region A. In order not to spoil the correct leading order 
for a quasi-zero-dimensional situation as given by $1,4 one has the additional requirement 
/(£ — » 1) ~ (1 — £) n with n > 2, which implies $3^ ~ 5 2n ~ 2 . But in region B £ varies between 
zero and one so that always $3^ ~ The function / must be nonnegative in this range, 
otherwise the functional would not be bounded from below [this happened in the original 
FMT QTOJ for which /(£) = (1/3 — £ 2 )/(87r)]. Then the argument of Sec. |III C| runs through 
and the asymptotic profiles will always be step functions. We conclude that an improved 
description of the high density solid is not possible within FMT if only the scalar and vector 
weighted densities are used, the tensor weight function of FMT2 is inevitable. On the other 
hand in FMT2 the behavior near close packing is exclusively determined by <p2 so that no 
conditions on the precise form of ^3 can be deduced. 



ACKNOWLEDGMENTS 

We thank Sander Pronk for providing the Monte Carlo program and Y. Rosenfeld for 
stimulating discussions. This work is part of the research program of the Stichting voor 
Fundamenteel Onderzoek der Materie (Foundation for Fundamental Research on Matter) 
and is made possible by financial support from the Nederlandse Organisatie voor Weten- 
schappelijk Onderzoek (Netherlands Organization for the Advancement of Research). B.G. 
acknowledges the financial support of the EU through the Marie Curie TMR Fellowship 
programme. 



APPENDIX A: GAUSSIAN PEAKS IN FMT1 



For Gaussian density peaks po(s) = it 3 / 2 e s2 the leading contributions to the weighted 
densities are (see Eqs. (f46|)-(|48|)) 

7i3o(t) = ~(1 - erf(t)) n 20 (t) = (Al) 

and for a width A = a/ (2a) the dominant excess free energy contribution is (see Eq. (|53|) ) 
62 a 



3B 



7T 



3/2 A 



dp / dz 



e (p+a+z) 2 + e ( p+a -z) 2 



-3 



X 



- erf (p + a + z) + - erf (p + a - z) 

— _ 



1 -2 



+ 0(A°). (A2) 
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In order that F ex does not become too large for a — > we expect a — > oo. In this limit the 
substitutions p' = pa and z' = za yield 



3B 



1 cr exp(— 3a 2 
12^ A ^ 



(A3) 



Now we can add 



^/iV = -^ln7r(A/A) 2 -^ 

and minimize with respect to A which gives 

cr 1 3 a 2 

A2^ eXp( -4A^ ) = 3 - 

This equation indeed has a solution with A/a —>■ for a — > 0; solved for a one has 



(A4) 



(A5) 



2A 

7! 



— In ( 6y/n 



-A 



1/2 



(A6) 



which demonstrates that A/a decays only very slowly. Nevertheless this decay is at variance 
with the physical expectation A/a — > const which is well supported by computer simulations 



TJL and, as shown in the main text, is also fulfilled within the present theory if allowance 



is made for more general profile shapes. 



17 



REFERENCES 



[1] S. Alexander and J. McTague, Phys. Rev. Lett. 41, 702 (1978). 
[2] B. Groh and B. Mulder, Phys. Rev. E 59, 5613 (1999). 
[3] Y. Singh, Phys. Rep. 207, 351 (1991). 

[4] R. Evans, in Fundamentals of Inhomogeneous Fluids, edited by D. Henderson (Marcel 

Dekker, New York, 1992), Chap. 3, pp. 85-175. 
[5] A. Haymet, in Fundamentals of Inhomogeneous Fluids, edited by D. Henderson (Marcel 

Dekker, New York, 1992), Chap. 9, pp. 363-405. 
[6] R. Ohnesorge, Ph.D. thesis, Universitat Miinchen, 1994. 

[7] N. Ashcroft, in Density Functional Theory, edited by E. Gross and R. Dreizler (Plenum, 

New York, 1995), pp. 581-623. 
[8] T. Ramakrishnan and M. Yussouff, Phys. Rev. B 19, 2775 (1979). 
[9] A. Haymet, J. Chem. Phys. 78, 4641 (1983). 
[10] Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989). 

[11] Y. Rosenfeld, M. Schmidt, H. Lowen, and P. Tarazona, Phys. Rev. E 55, 4245 (1997). 
[12] R. Ohnesorge, H. Lowen, and H. Wagner, Europhys. Lett. 22, 245 (1993). 
[13] B. Laird, J. McCoy, and A. Haymet, J. Chem. Phys. 87, 5449 (1987). 
[14] D. Young and B. Alder, J. Chem. Phys. 60, 1254 (1974). 

[15] J. P. Hansen and I. McDonald, Theory of simple liquids, 2nd ed. (Academic, London, 
1986). 

[16] M. Baus and J. Colot, Mol. Phys. 55, 653 (1985). 
[17] M. Baus, J. Phys. Condens. Matter 1, 3131 (1989). 
[18] Y. Rosenfeld, J. Phys. Condens. Matter 8, L795 (1996). 
[19] Y. Rosenfeld and P. Tarazona, Mol. Phys. 95, 141 (1998). 
[20] R. Kayser and H. Raveche, Phys. Rev. A 17, 2067 (1978). 
[21] J. Scheurle, J. Math. Anal. Appl. 59, 596 (1977). 

[22] R. Buehler, R. Wentorf, J. Hirschfelder, and C. Curtiss, J. Chem. Phys. 19, 61 (1950). 

[23] J. Kirkwood, J. Chem. Phys. 18, 380 (1950). 

[24] W. Wood, J. Chem. Phys. 20, 1334 (1952). 

[25] W. Hoover, J. Chem. Phys. 43, 371 (1965). 

[26] Z. Salsburg and W. Wood, J. Chem. Phys. 37, 798 (1962). 

[27] B. Alder, W. Hoover, and D. Young, J. Chem. Phys. 49, 3688 (1968). 

[28] P. Tarazona and Y. Rosenfeld, Phys. Rev. E 55, R4873 (1997). 

[29] C. Tejero, M. Ripoll, and A. Perez, Phys. Rev. E 52, 3632 (1995). 

[30] C. Tejero, Phys. Rev. E 55, 3720 (1997). 



18 



FIGURES 



FIG. 1. The left-hand side of Eq. (||) for the hard-sphere direct correlation function in the Per- 
cus-Yevick approximation. The wave number k max (p) corresponds to the maximum of cpy(p, k) at 
a given density p = p*a~ 3 . The curve lies below unity for all admissable densities p* < 6/ir = 1.910, 
i.e., for packing fractions rj = p*ir/6 < 1, which means that there is no liquid-solid spinodal. The 
close-packing limit occurs at p* = y2. 

FIG. 2. Density profiles in a high density fee crystal calculated from Ramakrishnan-Yussouff 
DFT. Note that the distance r from the lattice site and the density are scaled by the free distance 
a = R nn — a, which varies over 2.5 orders of magnitude in this density range. 

FIG. 3. The same as Fig. || but using p* = 0.95 as the density argument of the DCF. In this 
case the profiles are monotonia 

FIG. 4. Widths A corresponding to minima (lower branches) and saddle points (upper 
branches) of the Ramakrishnan-Yussouf functional restricted to Gaussian profiles for fee and bec 
solids. The asymptotic linear behavior indicated by the dashed lines was calculated from Eq. (|2~T|). 

FIG. 5. Free energies per particle of high density solids from the Ramakrishnan-Yussouff, the 
two versions of fundamental-measure DFT, and from molecular dynamics |27j. The de Broglie wave 
length has been set to the particle diameter. The asymptotic behavior indicated by the dotted lines 
is logarithmic in the free distance a (see Eq. (^2|)) in all cases. 

FIG. 6. The space between two nearest-neighbor sites (black dots) in a crystal. The radii of 
the spheres is a/2 ± O(A) where A is the width of the (spherical) density peaks. The weighted 
densities in region A are only influenced from one site while in B both sites contribute. In the 
remaining space the excess free energy densities (pi are neglible. 

FIG. 7. Density profiles obtained from the fundamental-measure theory FMT1. 

FIG. 8. Density profiles obtained from the improved fundamental-measure theory FMT2. 

FIG. 9. Illustration of the cells C for which the fundamental measure theory predicts a constant 
density (a 2d analogon of the 3d crystal is drawn). The circular arcs and their straight connection 
limit the set of points whose distance to C is smaller than a/2, i.e. the region A. The point P 
cannot belong to C because otherwise P' would belong to C and some points in between would 
have distances smaller than a/2 from both C and C, i.e. region B would be nonempty. 
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FIG. 10. Comparison of the density profiles from Monte-Carlo simulation and the den- 
sity-functional theory FMT2 for two bulk densities. A Gaussian profile would correspond to a 
straight line in this plot. 
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